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Abstract 


This  paper  evaluates  the  performance  of  various  orbit  propagation  theo¬ 
ries  for  artificial  earth  satellites  in  different  orbital  regimes.  Specifically,  R&D 
GTDS’s  Cowell  (numerical  technique),  DSST  (semianalytical  technique),  SGP, 
SGP4,  and  Brouwer-Lyddane  (analytic  techniques)  orbit  propagators  are  com¬ 
pared  for  decaying  circular  (~200  km  perigee  height),  low  altitude  circular  (590 
km  perigee  height),  high  altitude  circular  (1340  km  perigee  height),  Molniya, 
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and  geosynchronous  orbits.  All  test  cases  implement  a  one  orbital  period  dif¬ 
ferential  correction  fit  to  simulated  data  derived  from  a  Cowell  truth  trajectory. 
These  fits  are  followed  by  a  one  orbital  period  predict  with  the  DC  solve-for 
vector.  Trajectory  comparisons  are  made  with  the  Cowell  “truth”  trajectory 
over  both  the  fit  and  predict  spans.  Computation  time  and  RMS  errors  are 
used  as  comparison  metrics.  The  Unix-based  version  of  R&D  GTDS  (NPS  SUN 
Sparc  10)  is  the  test  platform  used  in  this  analysis. 


Nomenclature 


a  Semimajor  Axis 

AOG  DSST  Averaged  Orbit  Generator 

BL  Brouwer- Lydanne  Analytic  Propagator 

B*  SGP4  Drag  Parameter 

C'd  Drag  Coefficient 

Cr  SRP  Parameter 

DC  Differential  Correction 

DSST  Draper  Semianalytic  Satellite  Theory 

e  Eccentricity 

E  Eccentric  Anomaly 

f  True  Anomaly 

F  Eccentric  Longitude  (E  +  u>  +  SI) 

GEM  Goddard  Earth  Model 

GPS  Global  Positioning  System 

HST  Hubble  Space  Telescope 

i  Inclination 

ITER  Number  of  DC  Iterations 

JGM  Joint  Gravitational  Model 

kg  Kilograms 

km  Kilometers 

L  True  Longitude  (f  +  ui  +  SI) 

m  Meters 

M  Mean  Anomaly 

min  Minutes 

MS1S  Mass  Spectrometer  Incoherent  Scatter 

n,m  Degree,  Order  of  Gravity  Field 

h/2  SGP  Drag  Parameter 

Npq  BL  Drag  Parameter 

OD  Orbit  Determination 

opt  Optimized 

ORB1  GTDS  Output  Ephemeris  File 

osc  Osculating  Elements 

RhD  GTDS  Draper  Laboratory’s  Version  of  the  Goddard  Trajectory  Determination  System 
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RMS 

Root  Mean  Square 

SATCAT 

Satellite  Catalog 

SCN 

Space  Control  Network 

sec 

Seconds 

SGP 

Analytic  Theory  Based  on  the  Work  of  Kozai13 

SGP4 

Analytic  Theory  Based  on  the  Work  of  Lane  and  Cranford 

SPG 

DSST  Short  Periodic  Generator 

SRP 

Solar  Radiation  Pressure 

SSN 

Space  Surveillance  Network 

TOPEX 

Topographical  Ocean  Experiment  Satellite 

WGS 

World  Geodetic  Survey 

A 

Mean  Longitude  (M  +  ui  +  ft) 

9 

Greenwich  Hour  Angle 

UJ 

Argument  of  Perigee 

ft 

Longitude  of  Ascending  Node 

Introduction 

R&D  GTDS  is  Draper  Laboratory’s  research-based  orbit  determination  testbed1 . 
This  analysis  tool  evolved  from  its  R&D  counterpart  at  the  Goddard  Space 
Flight  Center.  Dr.  Paul  Cefola,  Program  Manager  at  Draper  Laboratory  and 
Lecturer  at  the  Massachusetts  Institute  of  Technology,  has  overseen  the  develop¬ 
ment  and  expansion  of  this  testbed  (by  a  team  of  scientists  at  Draper  Laboratory 
and  a  continuing  string  of  graduate  students  at  MIT)  over  the  past  twenty  years 
(see  Fonte7).  In  its  current  form,  R&D  GTDS  is  capable  of  performing 

•  early  orbit  determination 

—  Gauss, 

—  Double-R  iteration, 

—  range  and  angles  methods, 

—  range-only  method, 

•  orbit  propagation 

—  numerical, 

—  analytical, 

—  seminanalytical  techniques, 

•  estimation 

—  batch  weighted  least  squares, 

—  filtering  algorithms, 
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•  error  analysis 


—  error  categorization  for  the  OD  process, 

•  data  simulation 

—  simulated  tracking  data  from  specified  ground  stations. 

This  paper  will  present  a  comparative  study  of  the  various  orbit  propaga¬ 
tion  techniques  (Cowell,  DSST,  SGP,  SGP4,  and  BL)  imbedded  in  R&D  GTDS. 
These  orbit  propagators  represent  a  standard  blend  of  special  and  general  per¬ 
turbation  techniques.  With  special  perturbation  techniques  (Cowell),  a  direct 
numerical  integration  of  the  equations  of  motion  is  implemented.  This  entails 
multiple  force  evaluations  to  “step”  from  a  given  set  of  initial  conditions  to  a 
final  solution.  If  the  initial  conditions  are  changed  or  altered,  new  evaluations 
must  be  made  at  each  of  the  various  time  steps.  Therefore,  each  solution  is 
unique  to  a  given  set  of  initial  conditions.  In  special  perturbation  methods, 
consideration  must  be  given  to  the  computation  time,  physical  accuracy,  trun¬ 
cation  error,  and  round-off  error  associated  with  selection  of  an  appropriate 
time  step  for  integration.  In  general,  high  accuracy  is  obtained  at  the  cost  of 
computational  efficiency  (these  techniques  have  typically  been  limited  to  “spe¬ 
cial  case”  processing  in  Space  Command’s  operational  environment  due  to  their 
inherent  computational  burden;  this  includes,  but  is  not  limited  to,  scenarios 
requiring  high  accuracy  or  a  limited  number  of  satellites) . 

General  perturbation  techniques  (SGP,  SGP4,  BL),  on  the  contrary,  do  not 
use  multiple  time  steps  to  transfer  from  a  set  of  initial  conditions  to  a  final 
solution.  Rather,  general  perturbation  techniques  implement  an  analytical  in¬ 
tegration  of  the  equations  of  motion  to  transform  a  set  of  initial  conditions 
directly  to  a  final  solution.  However,  due  to  the  complex  nature  of  the  equa¬ 
tions  representing  the  physical  models,  exactly  integrable  expressions  are  diffi¬ 
cult  to  obtain.  For  this  reason,  simplifications,  approximations,  and  truncations 
are  made  to  the  equations  of  motion  to  obtain  expressions  which  are  integrable. 
These  simplifications,  however,  greatly  reduce  the  accuracy  attainable  with  gen¬ 
eral  perturbation  techniques.  In  summary,  computational  efficiency  is  gained 
at  the  cost  of  accuracy  (these  techniques  have  typically  been  used  for  the  bulk 
of  Space  Command’s  operational  missions  due  to  the  computational  efficiency 
advantages). 

A  more  recent  approach  to  perturbation  analysis  is  semianalytic  techniques. 
These  techniques  combine  the  primary  advantageous  aspect  of  special  pertur¬ 
bation  methods  (high  accuracy)  with  that  of  general  perturbation  techniques 
(computational  efficiency).  For  one  such  semianalytic  technique,  the  Draper 
Semianalytic  Satellite  Theory  (DSST),  the  equations  of  motion  are  separated 
into  two  distinct  categories  (see  McClain3,4,  Danielson  et  al2).  One  category 
contains  the  secular  and  long-period  perturbative  contributions  to  a  satellite’s 
motion;  the  portion  of  the  software  which  propagates  this  motion  is  referred  to  as 
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the  averaged  orbit  generator  (AOG),  and  is  based  upon  “mean”  elements.  The 
other  category  contains  the  short-period  perturbative  contributions  to  a  satel¬ 
lite’s  motion,  which  are  modeled  in  the  short  periodic  generator  (SPG;  DSST 
represents  the  state  of  the  art  in  short  periodic  models).  The  separation  of  the 
short-periodic  contributions  from  the  secular  and  long-period  contributions  is 
accomplished  via  the  generalized  method  of  averaging  (see  Morrison9).  With 
the  short-periodic  contributions  separated  (which,  due  to  their  high  frequency 
nature,  are  step-size  constraining),  the  averaged  equations  of  motion  can  be 
integrated  numerically  with  large  step  sizes  (typically  on  the  order  of  a  day), 
while  the  short  periodics  can  be  recovered  analytically  at  output  times  (with  the 
use  of  Fourier  analysis  and  interpolation  schemes).  The  specific  force  models 
included  in  the  AOG  and  SPG  can  be  chosen  by  the  user  at  run  time,  creating  a 
theory  that  is  highly  accurate,  efficient,  and  flexible.  It  is  of  the  authors’  opinion 
that  these  techniques  are  optimally  suited  for  operational  processing  (as  well  as 
for  a  multitude  of  other  applications,  such  as  long-term  mission  planning). 

With  special,  general,  and  semianalytic  perturbation  techniques  defined,  it 
is  now  possible  to  describe  the  specific  force  models  available  in  each  of  the 
propagators  analyzed  in  this  study: 

Cowell  *  (user  has  ability  to  select  force  models)  (see  GTDS  Mathematical 
Specification1) 

•  Geopotential 

—  models  ranging  from  GEM1  to  GEMT3,  JGM1,  JGM2,  WGS72,  and 
WGS84 

—  maximum  degree  and  order  of  50  (n,m  <  50)  t 

•  Atmospheric  Drag 

—  Harris-Priester  (static  atmosphere;  tabulated  values  of  density  versus 
altitude) 

—  Jacchia- Roberts  (dynamic  atmosphere;  solar  flux  and  geomagnetic 
indices) 

—  Jacchia  70 

—  Jacchia  64 

—  MSIS 


*It  should  be  noted  that  the  same  force  models  are  available  for  both  Cowell  and  DSST. 
Geopotential  models,  atmosheric  data,  solar  lunar  and  planetary  ephemerides,  timing  coeffi¬ 
cients,  and  Newcomb  operator  data  exist  in  a  binary  data  base  attached  to  R&D  GTDS.  As 
new  models  or  raw  inputs  become  available,  the  binary  data  bases  are  updated  and  easily 
attached  to  R&D  GTDS. 

^50  by  50  software  is  available  in  a  library  linkable  to  R&D  GTDS.  The  current  code  at 
Naval  Postgraduate  School  has  a  maximum  degree  and  order  of  21 
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•  Third-Body  (Lunar-Solar)  Point  Mass  (JPL  solar,  lunar,  and  planetary 
ephemerides;  extendable  to  multi-body) 

•  Solar  Radiation  Pressure  (cylindrical  shadow  model) 


DSST  +  (user  has  ability  to  select  force  models)  (see  McClain3,4,  and  Daniel¬ 
son  et  al  2) 

•  Geopotential 

—  models  ranging  from  GEM1  to  GEMT3,  JGM1,  JGM2,  WGS72  and 
WGS84 

—  maximum  degree  and  order  of  50  (n,m  <  50)  § 

—  First  Order,  AOG 

*  recursive,  closed  form  zonal  model 

*  recursive  tesseral  resonance  model 
—  First  Order,  SPG 

*  recursive,  closed  form  zonal  short  periodic  model  (Fourier  series 
in  L) 

*  recursive,  closed  form  tesseral  m-daily  short  periodic  model  (Fourier 
series  in  6) 

*  recursive  tesseral  linear  combination  short  periodic  model  (Fourier 
series  in  A  and  6) 

—  Second  Order,  AOG 

*  </|  to  first  order  in  e 
—  Second  Order,  SPG 

*  </|  to  zeroth  order  in  e 

—  Second  Order  Coupling  Terms,  AOG 

*  /  drag  coupling 

—  Second  Order  Coupling  Terms,  SPG 

*  recursive,  closed  form  J 2  /  tesseral  m-daily  model 

*  Drag  /  m-daily 


+  It  should  be  noted  that  the  same  force  models  are  available  for  both  Cowell  and  DSST. 
Geopotential  models,  atmosheric  data,  solar  lunar  and  planetary  ephemerides,  timing  coeffi¬ 
cients,  and  Newcomb  operator  data  exist  in  a  binary  data  base  attached  to  R&D  GTDS.  As 
new  models  or  raw  inputs  become  available,  the  binary  data  bases  are  updated  and  easily 
attached  to  R&D  GTDS. 

§50  by  50  software  is  available  in  a  library  linkable  to  R&D  GTDS.  The  current  code  at 
Naval  Postgraduate  School  has  a  maximum  degree  and  order  of  21 
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•  Atmospheric  Drag,  AOG  and  SPG 

—  same  drag  models  as  described  for  Cowell 

—  numerical  averaging  (quadrature) 

—  SPG  Fourier  series  expressed  in  A 

—  some  optional  second  order  drag  effects 

•  Third-Body  Point  Mass  (JPL  solar,  lunar,  and  planetary  ephemerides; 
extendable  to  multi-body) 

—  first  order  AOG 

*  recursive,  closed  form  lunar-solar  model 

—  first  order  SPG 

*  recursive,  closed  form  lunar-solar  model  (Fourier  series  in  F) 

*  weak  time  dependence  corrections 

—  double  averaged  formulation 

•  Solar  Radiation  Pressure,  AOG  and  SPG 

—  cylindrical  shadow  model 

—  numerical  averaging 

—  SPG  Fourier  series  expressed  in  A 


SGP  (no  force  model  flexibility  other  than  drag  parameter  adjustment)  (see 
Floots5,  Fieger6,  and  Flerriges10) 

•  gravitational  model  based  on  a  simplification  of  the  work  of  Kozai13 

—  J2  secular  effects 
—  J2  and  J3  long-period  effects 

—  truncated  J2  short  periodics  to  zeroth  order  in  e 

—  secular  drag  effects  in  a,  e,  and  M  (drag  model  based  on  Taylor  series 
in  mean  anomaly  (h/ 2);  drag  effect  modeled  on  eccentricity  such  that 
perigee  height  remains  constant). 


SGP4  (no  force  model  flexibility  other  than  drag  parameter  adjustment)  (see 
Floots5,  Fieger6,  Flerriges10,  Lane  and  Floots11,  and  Flujsak  and  Floots12) 

•  gravitational  model  based  on  a  simplification  of  the  work  of  Brouwer14 

—  J2,  J4,  and  j|  secular  effects 
—  J2  and  J3  long-period  effects 
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—  truncated  J2  short  periodics  to  zeroth  order  in  e 

—  tesseral  terms  (2,2),  (3,2),  (5,2),  (4,4),  and  (5,4)  for  orbits  with  peri¬ 
ods  >  225  minutes;  designed  for  geosynchronous  and  Molniya  orbits 

—  first  order  lunar-solar  point  mass  terms  for  orbits  with  periods  >  225 
minutes 

—  drag  model  based  on  power  density  function  ((q0  —  s)4,  B*),  secular 
effects  only 


R&D  GTDS  BL  (no  force  model  flexibility)  (see  GTDS  Mathematical  Specification1) 

•  Lyddane  modified  Brouwer  theory  to  obtain  algorithms  applicable  for  sin¬ 
gularities  (i  =  0  or  e  =  0) 

—  J 2 ,  J4,  and  j|  secular  effects 
—  J2 ,  J3 ,  J4,  and  J5  long-period  effects 
—  J2  short  periodics 

—  drag  model  based  on  Brouwer  drag  coefficients,  Np  q  (“free  parame¬ 
ters”  which  can  be  solved-for  in  the  DC  to  better  model  drag);  correct 
mean  anomaly  only 


R&D  GTDS  OD  Options 

•  Coordinate  frames 

—  mean  Earth  equator  and  equinox  of  1950.0 
—  FK4  derived  true  of  date  (OhOmOs  zulu) 

—  mean  Earth  equator  and  equinox  of  2000.0  (J2000;  available  through 
externally  interfaced  program) 

—  FK5  derived  true  of  date  (true  of  time  of  interest;  available  through 
externally  interfaced  program) 

•  Time  standards 

—  implements  UTC  (broadcast  time)  as  standard  time 

—  other  time  standards  available  via  conversion  utilities 
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Test  Protocol 


The  comparisons  in  this  analysis  are  based  against  Cowell  truth  trajectories; 
one  Cowell  truth  trajectory  was  established  for  each  orbital  regime.  For  the 
low  altitude  circular,  high  altitude  circular,  Molniya  and  geosynchronous  test 
cases,  the  osculating  elements  (and  appropriate  perturbation  parameters)  used 
to  generate  the  truth  trajectories  were  derived  from  Cowell  fits  to  an  SGP4 
ephemeris  derived  from  actual  two-card  element  sets  (SATCAT).  For  the  decay 
case,  the  osculating  elements  were  taken  from  Dyar17.  This  specific  information 
pertinent  to  the  truth  trajectories  is  as  follows: 


Epoch  Date 

Epoch  Time  (hhmmss.ssss) 
a  (osc) 
e  (osc) 
i  (osc) 

SI  (osc) 
u>  (osc) 

M  (osc) 

Period 

Perigee  Eleight 
Gravity  Model 
Atmos.  Dens. 

Third  Body 
SRP 

Integrator 
Step  Size 

Input,  Output  and  Integration  Frame 
Spacecraft  Area,  Mass 
ORB1  Output  Frequency 


24  February  82 
00:00:00.0000 
6628.45287  km 
0.008921 
64.83828° 

224.50715° 

271.80047° 

115.250818° 

89.5113  min 
191.1834  km 
21x21  GEM10B 

Jacchia- Roberts  Drag  (Cd  =  2.0) 

Lunar  -  Solar  Point  Mass 
No 

12th  Order  Summed  Cowell/  Adams  Predict-Partially  Correct 
60  sec 

Mean  Earth  Equator  and  Equinox  of  1950.0 
1  m2,  100  kg 
Every  450  sec 


Table  1:  Decaying  Orbit 
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Epoch  Date 

Epoch  Time  (hhmmss.ssss) 
a  (osc) 
e  (osc) 
i  (osc) 

SI  (osc) 
u>  (osc) 

M  (osc) 

Period 

Perigee  Height 
Gravity  Model 
Atmos.  Dens. 

Third  Body 
SRP 

Integrator 
Step  Size 

Input,  Output  and  Integration  Frame 
Spacecraft  Area,  Mass 
0RB1  Output  Frequency 


4  December  94 

21:10:59.5456 

6980.108566 

0.00170791 

28.327265° 

218.492078° 

40.611761° 

319.915316° 

96.7270  min 
590.05  km 
21x21  GEM10B 

Jacchia- Roberts  (Schatten)  Drag  (Cd  =  2.59856) 

Lunar  -  Solar  Point  Mass 
CR  =  8.68846 

\2th  Order  Summed  Cowell/  Adams  Predict-Partially  Correct 
60  sec 

Mean  Earth  Equator  and  Equinox  of  1950.0 
1  m2,  100  kg 
Every  400  sec 


Table  2:  Low  Altitude,  Circular  Orbit  (SATCAT  #  20580,  HST) 


Epoch  Date 

Epoch  Time  (hhmmss.ssss) 
a  (osc) 
e  (osc) 
i  (osc) 

0  (osc) 
u>  (osc) 

M  (osc) 

Period 

Perigee  Height 
Gravity  Model 
Atmos.  Dens. 

Third  Body 
SRP 

Integrator 
Step  Size 

Input,  Output  and  Integration  Frame 
Spacecraft  Area,  Mass 
ORB1  Output  Frequency 


4  December  94 

11:30:01.5434 

7721.500966 

0.000482180 

65.830824° 

296.770687° 

19.611311° 

340.277743° 

112.5425  min 
1339.6408  km 
21x21  GEM10B 

Jacchia- Roberts  (Schatten)  Drag  (C'd  =  1.8661152) 

Lunar  -  Solar  Point  Mass 
CR  =  -5.4211956 

12th  Order  Summed  Cowell/  Adams  Predict-Partially  Correct 
60  sec 

Mean  Earth  Equator  and  Equinox  of  1950.0 
1  m2,  100  kg 
Every  460  sec 


Table  3:  High  Altitude,  Circular  Orbit  (SATCAT  #  22076,  TOPEX) 
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Epoch  Date 

Epoch  Time  (hhmmss.ssss) 
a  (osc) 
e  (osc) 
i  (osc) 

SI  (osc) 
u>  (osc) 

M  (osc) 

Period 

Perigee  Height 
Gravity  Model 
Atmos.  Dens. 

Third  Body 
SRP 

Integrator 
Step  Size 

Input,  Output  and  Integration  Frame 
Spacecraft  Area,  Mass 
0RB1  Output  Frequency 


29  August  94 

06:25:23.8964 

26586.72673 

0.731526 

63.958664° 

47.836839° 

289.061990° 

9.509759° 

~12  hours 
759.70787  km 
21x21  GEM10B 

Jacchia- Roberts  (Schatten)  Drag  (Cr  =  160.9959772) 

Lunar  -  Solar  Point  Mass 
CR  =  4.5991303 

\2th  Order  Summed  Cowell/  Adams  Predict-Partially  Correct 
60  sec 

Mean  Earth  Equator  and  Equinox  of  1950.0 
1  m2,  100  kg 
Every  450  sec 


Table  4:  Molniya  Orbit  (SATCAT  #  21897) 


Epoch  Date 

Epoch  Time  (hhmmss.ssss) 
a  (osc) 
e  (osc) 
i  (osc) 

0  (osc) 
u>  (osc) 

M  (osc) 

Period 

Perigee  Height 
Gravity  Model 
Atmos.  Dens. 

Third  Body 
SRP 

Integrator 
Step  Size 

Input,  Output  and  Integration  Frame 
Spacecraft  Area,  Mass 
ORB1  Output  Frequency 


31  August  94 

16:19:28.4739 

42549.449791 

0.000487508 

11.162079° 

47.617617° 

56.938285° 

74.862562° 

~24  hours 
36150.5696  km 
21x21  GEM10B 
No 

Lunar  -  Solar  Point  Mass 
CR  =  5.2927536 

12th  Order  Summed  Cowell/  Adams  Predict-Partially  Correct 
60  sec 

Mean  Earth  Equator  and  Equinox  of  1950.0 
1  m2,  100  kg 
Every  450  sec 


Table  5:  Geosynchronous  Orbit  (SATCAT  #  10516) 
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Note:  Even  though  the  orbits  used  in  this  analysis  can  be  tied  to  a  specific 
object  number  in  the  SATCAT,  a  standard  area  of  1  m2  and  mass  of  100  kg 
were  implemented  throughout  the  test  cases  (which  provides  a  “rule-of-thumb” 
area  to  mass  ratio).  In  addition,  all  test  cases  except  the  decaying  orbit  have 
epochs  in  1994;  the  decay  case,  which  was  chosen  to  demonstrate  the  specific 
effects  of  atmospheric  drag,  used  a  “noisy”  epoch  in  1982  (geomagnetic  indices 
were  disturbed  in  this  time  period). 

The  Jacchia- Roberts  (Schatten)  hie  represents  a  particular  version  of  a  Jacchia- 
Roberts  hie  generated  by  Ken  Schatten  (Goddard  Space  Flight  Center)  and 
David  Carter18  (Draper  Laboratory).  Schatten  uses  monthly  values  of  solar 
hux  and  geomagnetic  activity  to  generate  a  “smooth”  (yet  still  dynamic)  atmo¬ 
sphere.  Carter  applied  interpolation  techniques  to  reduce  Schatten’s  monthly 
values  to  “smooth”  daily  values,  which  are  required  by  R&D  GTDS.  Schatten’s 
hies  performed  as  well  as  any  other  atmospheric  model  during  the  last  peak  of 
the  solar  cycle,  and  for  this  reason,  have  been  included  in  the  binary  data  base 
of  R&D  GTDS  (Cefola19).  More  details  concerning  these  procedures  can  be 
found  in  the  work  of  Sabol15. 

All  tests  in  support  of  this  analysis  implemented  one-orbital  period  DC  hts 
to  the  appropriate  Cowell  truth  model.  The  “observations”  used  in  these  hts 
were  simulated  data  in  the  form  of  GTDS  ORB1  hies  (evenly  spaced,  time 
tagged  values  of  position  and  velocity;  very  much  like  GPS  navigation  solu¬ 
tions)  based  on  the  Cowell  trajectories.  One  orbital  period  predictions  with 
the  DC  solve-for  vectors  were  then  generated  and  compared  with  the  Cowell 
truth  trajectories.  RMS  errors  over  the  ht  (DC  RMS)  and  predict  (R&D  GTDS 
Ephemeris  Comparison  Program  RMS)  spans  were  used  as  accuracy  metrics; 
timing  comparisons  were  based  on  a  call  to  an  internal  clock  routine  at  the 
initiation  and  termination  of  the  program. 

Results  of  Testing 


For  each  of  the  hve  classes  of  satellite  orbits  (decaying,  low  altitude  circu¬ 
lar,  high  altitude  circular,  Molniya  and  geosynchronous)  several  orbit  propaga¬ 
tors/orbit  propagator  conhgurations  were  analyzed.  In  general,  the  following 
protocol  was  used: 

•  One  Cowell  DC  with  force  models  matching  those  in  the  truth  ephemeris; 
this  run  ensures  the  DC  can  reconstruct  the  truth  ephemeris  from  per¬ 
turbed  initial  conditions  (i.e. ,  test  if  the  DC  is  functioning  properly). 

•  Cowell  DCs  with  force  models  truncated  as  compared  to  the  truth  ephemeris. 
These  runs  analyze  the  “pure”  impact  of  truncating  force  models  (i.e.,  re¬ 
ducing  the  geopotential  from  21  x  21  to  4  x  4,  etc.). 

•  One  DSST  DC  with  force  models  configured  to  balance  computational 
efficiency  and  accuracy;  these  conhgurations  are  based  on  the  work  of 
Fonte  and  Sabol16. 
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•  Various  analytic  DCs  (BL,  SGP  and  SGP4)  with  and  without  drag  solve- 
fors.  It  should  be  noted  solving  for  drag  parameters  in  orbital  regimes 
where  drag  doesn’t  have  a  large  effect  was  intentional;  in  these  cases, 
the  drag  parameter  simply  represents  a  “free-parameter”  in  the  DC  at¬ 
tempting  to  absorb  errors  stemming  from  truncations  in  force  models  for 
perturbations  other  than  drag. 

In  addition,  several  different  element  sets  were  used  as  initial  conditions  for 
the  DC  in  each  case.  With  the  exception  of  the  decay  case,  results  quoted  in 
this  paper  are  based  on  using  the  two-card  element  set  from  which  the  truth 
trajectories  were  generated.  This  methodolgy  provided  a  standard  set  of  initial 
conditions  from  which  to  analyze  results.  It  should  be  noted  this  may  provide 
an  unfair  advantage  to  SGP  and  SGP4  (i.e. ,  Cowell  uses  osculating  elements  and 
therefore,  the  DC  may  have  to  execute  extra  iterations  in  order  to  converge). 
This  protocol  must  be  considered  when  analyzing  timing  results  (the  initial 
conditions  for  the  decaying  orbit  DCs  were  taken  from  osculating  truth  output 
ephemeris). 

Furthermore,  three  sets  of  timing  statistics  given  for  each  case:  time  to  per¬ 
form  the  DC  only,  time  to  perform  the  DC  and  subsequent  ephemeris  generation, 
and  DC  time  per  iteration.  The  total  time  to  perform  the  DC  is  representative 
of  processing  times  for  operational  catalog  maintenance;  ephemeris  generation 
times  provide  insight  to  the  speed  of  pure  orbit  propagation  with  the  various 
theories.  DC  time  per  iteration  provides  an  alternative  comparison  metric  for 
theories  which  require  a  different  number  of  iterations  to  arrive  at  a  final  so¬ 
lution.  However,  processing  is  based  upon  arriving  at  a  final  solution,  so  both 
time  per  DC  iteration  and  total  DC  time  are  meaningful  metrics. 

Results  can  be  now  given  for  each  of  the  various  orbit  classes. 

DECAYING  ORBIT: 

Configuration  Notes:  In  this  case,  Cowel^i  has  force  models  matching  those 
in  the  truth  (21x21  GEM10B  geopotential,  Jacchia- Roberts  drag,  lunar  solar 
point  mass  effects,  solve-for  drag  parameter  in  DC).  CowelU  is  identical  to 
Cowell2i,  with  the  exception  of  a  reduced  geopotential  (i.e.,  4x4).  Cowell;) 
differs  from  Cowell  only  in  that  the  exact  drag  term  from  the  truth  ephemeris 
has  been  hardwired  into  the  run  (it  is  not  being  solved  for  by  the  DC).  The 
annotation  of  Npq  or  “drag”  after  a  theory  indicates  DC  solve-for  parameters 
used  in  the  run. 

The  results  given  in  Table  6  clearly  show  the  dominance  of  atmospheric  and 
gravitational  effects.  For  example,  reducing  the  Cowell  geopotential  configura¬ 
tion  from  21x21  to  4x4  results  in  ~265  m  predict  error  over  one  revolution  (refer 
to  Cowell;)).  The  CowelU  results  highlight  the  DC  “observability”  problems  as¬ 
sociated  to  coupling  an  extremely  short  fit  span  with  theories  having  severely 
truncated  force  models;  the  CowelU  DC,  which  attempted  to  solve  for  a  drag 
term,  produced  an  error  ~100  m  larger  than  Cowell;)  (~375  m).  Inspection  of 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 

RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

BL 

229.40839 

6 

5432.4 

1.88 

1.37 

.23 

BL  (Npq) 

236.62321 

5 

2094.90 

2.75 

1.32 

.26 

CowelCi 

2.60343(E-9) 

5 

2.2922(E-8) 

9.17 

2.77 

.55 

Cowell 

22.761151 

4 

374.11 

6.96 

2.03 

.51 

Cowell;) 

22.292897 

4 

264.14 

7.74 

2.34 

.59 

SGP4 

138.36233 

4 

8897.6 

1.96 

1.45 

.36 

SGP4  (drag) 

138.63008 

5 

503.23 

2.56 

1.66 

.33 

SGP 

638.31315 

4 

7858.4 

2.46 

1.40 

.35 

SGP  (drag) 

1182.5092 

4 

12598.0 

2.47 

1.14 

.29 

DSST  (opt) 

17.416619 

6 

65.497 

6.83 

5.17 

.86 

Table  6:  Decay  Orbit  90  Minutes  Fit  and  90  Minutes  Predict 


the  solved-for  drag  term  in  the  CowelU  DC  indicated  the  DC  could  not  exactly 
recover  the  drag  term  used  in  the  generation  of  the  truth  ephemeris  (remember, 
the  Cowell;)  DC  used  the  exact  drag  term  used  in  the  generation  of  the  truth 
orbit  in  order  to  determine  the  “pure”  geopotential  error).  The  analytic  theories 
also  experience  these  observability  problems,  which  are  further  compounded  by 
their  lack  of  tesseral  terms  and  simplistic  drag  models.  Operationally,  these 
problems  are  countered  by  increased  fit  spans  in  an  attempt  to  “observe”  more 
of  the  effect  which  has  been  neglected  (for  decay  cases,  in  which  it  may  not  be 
practical  or  possible  to  increase  the  fit  span;  increased  amounts  of  observational 
data  can  also  be  used).  In  this  manner,  the  solve- for  parameters  in  the  DC  ab¬ 
sorb  some  of  the  error  stemming  from  truncated  force  models  (the  parameters, 
therefore,  become  less  physical;  even  though  the  fit  may  be  improved,  predic¬ 
tions  with  DC  solve- for  vector  become  worse).  The  optimized  DSST  results 
indicate  properly  configured  perturbation  theories  can  still  produce  excellent 
results  with  the  short  fit  spans;  however,  the  timing  statistics  suggest  that 
further  refinement  of  the  optimized  DSST  could  decrease  the  amount  of  time 
required  for  the  run-but  at  the  expense  of  accuracy  (users  can  configure  DSST 
for  specific  requirements). 

It  is  worth  stressing  atmospheric  effects  are  the  limiting  factor  in  the  decay 
case.  Further  testing  with  Cowell  and  DSST  was  undertaken  in  an  attempt 
to  more  realistically  represent  real  world  atmospheric  conditions.  These  DCs 
used  a  static  drag  model  (Harris-Priester)  to  fit  the  dynamic  (Jacchia- Roberts) 
truth  model.  This  protocol  takes  into  account  that  Jacchia- Roberts  atmospheric 
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parameters  (F10.7  solar  flux  values  and  geomagnetic  indices)  are  very  difficult 
to  predict  or  may  not  be  obtained  in  a  timely  manner.  Therefore,  fitting  the 
noisy  Jacchia- Roberts  model  with  the  smooth  Harris-Priester  model  (hottest 
H-P  table)  emulates  a  lack  of  perfect  knowledge  of  the  atmospheric  parameters. 
The  results  of  this  testing  are  given  in  Table  7. 


Theory 

DC 

RMS  (m) 

ITER 

Predict 
RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

Cowell2i 

1.190631 

4 

183.30 

7.53 

2.06 

.52 

DSST  (opt) 

16.308181 

6 

149.86 

5.07 

3.74 

.62 

Table  7:  Decaying  Orbit  (Real  World)  90  Minutes  Fit  and  90  Minutes  Predict 

These  results  indicate  the  optimized  DSST  performs  as  well  as  Cowell2i  with 
a  significant  savings  in  computational  time. 

LOW  ALTITUDE  CIRCULAR  ORBIT: 

Configuration  Notes:  In  this  case,  Cowell2i  has  force  models  matching  those 
in  the  truth  (21x21  GEM10B  geopotential,  Jacchia- Roberts  (Schatten)  drag, 
lunar  solar  point  mass  effects,  SRP,  solve-for  drag  and  SRP  parameters  in  DC). 
Cowell  is  identical  to  Cowell  i,  with  the  exception  of  a  reduced  geopotential 
(i.e. ,  4x4).  Cowell;)  differs  from  CowelU  only  in  that  the  exact  drag  and  SRP 
terms  from  the  truth  ephemeris  have  been  hardwired  into  the  run  (they  are  not 
being  solved  for  by  the  DC).  The  annotation  of  Npq  or  “drag”  after  a  theory 
indicates  DC  solve-for  parameters  used  in  the  run. 

In  this  analysis,  reducing  the  Cowell  geopotential  configuration  from  21x21 
to  4x4  results  in  ~42  m  predict  error  over  one  revolution  (refer  to  Cowell;)). 
Again,  DC  “observability”  problems  associated  to  coupling  a  short  fit  span  with 
theories  having  severely  truncated  force  models  is  evident  in  the  CowelU  DC; 
this  DC,  which  attempted  to  solve  for  both  a  drag  and  SRP  term,  produced  an 
error  ~160  m  larger  than  Cowell;)  (~200  m).  Inspection  of  the  solved-for  drag 
and  SRP  terms  in  the  Cowell  DC  indicated  the  DC  could  not  exactly  recover 
the  original  terms  used  in  the  generation  of  the  truth  ephemeris  (remember,  the 
Cowell;)  DC  used  the  exact  drag  and  SRP  terms  used  in  the  generation  of  the 
truth  orbit  in  order  to  determine  the  “pure”  geopotential  error).  In  addition, 
this  orbit  is  at  an  altitude  for  which  SRP  has  a  very  minimal  effect,  which  also 
limits  the  ability  of  the  DC  to  properly  “observe”  the  perturbing  effect  (in  fact, 
additional  testing  indicates  a  solution  for  the  SRP  term  is  not  required  at  this 
altitude  except  for  scenarios  requiring  the  utmost  accuracy) .  As  with  the  decay 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 

RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

BL 

50.287918 

4 

197.18 

1.84 

1.10 

.27 

BL  (Npq) 

49.534240 

4 

231.34 

1.88 

1.06 

.26 

Cowell2i 

2.72860(E-8) 

5 

0.19385(E-3) 

4.25 

3.00 

.60 

Cowell 

14.536095 

5 

202.50 

3.77 

2.49 

.50 

Cowell;) 

15.768104 

5 

41.655 

3.68 

2.88 

.58 

SGP4 

54.361123 

3 

196.34 

2.15 

1.03 

.34 

SGP4  (drag) 

53.842927 

4 

669.41 

2.35 

1.57 

.39 

SGP 

210.20938 

3 

933.95 

1.95 

1.06 

.35 

SGP  (drag) 

688.68465 

4 

14701.0 

2.07 

1.16 

.29 

DSST  (opt) 

34.889745 

4 

47.274 

2.99 

2.03 

.51 

Table  8:  HST  Orbit  100  Minutes  Fit  and  100  Minutes  Predict 


Theory 

DC 

RMS  (m) 

ITER 

Predict 
RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

CowelGi 

7.12470788(E-4) 

10 

0.35502 

5.60 

4.87 

.49 

Cowell 

14.595436 

5 

242.54 

3.75 

2.67 

.53 

DSST  (opt) 

34.862539 

4 

48.372 

2.71 

1.68 

.42 

Table  9:  HST  Orbit  (Real  World)  100  Minutes  Fit  and  100  Minutes  Predict 
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case,  the  performance  of  the  analytic  theories  is  limited  by  the  truncation  of 
the  geopotential  and  simplistic  atmospheric  models;  however,  this  orbit  class 
(as  well  as  the  high  altitude  circular  case)  represents  an  orbital  regime  in  which 
the  best  performance  from  the  analytic  theories  can  be  expected  .  Furthermore, 
observability  problems  can  also  be  noticed  with  the  analytic  theories  in  that 
solving  for  drag  terms  in  the  DC  produces  larger  errors  than  for  cases  without 
drag  parameter  solutions.  The  optimized  DSST  DC  again  performed  very  well 
in  terms  of  accuracy,  but  run  times  were  not  as  optimistic  as  expected.  To 
more  realistically  assess  the  run  time  for  optimized  DSST,  and  also  to  simulate 
operational  atmospheric  difficulties,  further  testing  was  accomplished  to  fit  the 
Flarris-Priester  model  (default  Fl-P  table)  to  the  Jacchia- Roberts  (Cowell)  truth. 
These  results  are  given  in  Table  9. 

These  results  more  accurately  reflect  the  blend  of  efficiency  and  accuracy 
that  can  be  obtained  with  DSST.  Further  results  for  this  orbit  class  can  be 
found  in  Fonte8,23. 

HIGH  ALTITUDE  CIRCULAR  ORBIT: 

Configuration  Notes:  In  this  case,  Cowel^i  has  force  models  matching  those 
in  the  truth  (21x21  GEM10B  geopotential,  Jacchia- Roberts  (Schatten)  drag, 
lunar  solar  point  mass  effects,  SRP,  solve-for  drag  and  SRP  parameters  in  DC). 
Cowell^  has  21x21  GEM10B  geopotential  and  lunar  solar  point  mass  effects 
only.  Cowell  is  identical  to  Cowell2i,  with  the  exception  of  a  reduced  geopo¬ 
tential  (i.e. ,  4x4).  Cowell;)  differs  from  CowelU  only  in  that  the  exact  drag  and 
SRP  terms  from  the  truth  ephemeris  have  been  hardwired  into  the  run  (they 
are  not  being  solved  for  by  the  DC).  The  annotation  of  Npq  or  “drag”  after  a 
theory  indicates  DC  solve-for  parameters  used  in  the  run. 

For  the  orbit  chosen  to  fulfill  the  high  altitude  circular  test  case,  reducing 
the  Cowell  geopotential  configuration  from  21x21  to  4x4  results  in  ~105  m  pre¬ 
dict  error  over  one  revolution  (refer  to  Cowell;)).  It  may  seem  peculiar  the  error 
using  a  4x4  geopotential  model  is  larger  for  this  orbit  than  for  the  lower  alti¬ 
tude,  HST  orbit.  This  discrepancy  can  be  attributed  to  the  larger  inclination 
for  the  TOPEX  orbit  (subsequent  testing  with  HST  at  an  inclination  similar  to 
the  TOPEX  orbit  produces  a  predict  error  ~  170m.)  As  with  the  decay  and 
low  altitude  test  cases,  DC  “observability”  problems  are  once  again  present; 
specifically,  an  error  ~150  m  larger  than  Cowell;)  is  obtained  when  the  DC  at¬ 
tempts  to  solve  for  drag  and  SRP  terms  in  the  Cowell  case  (predict  error  ~250 
m).  In  addition,  this  orbit  is  at  an  altitude  for  which  drag  has  a  small  effect, 
which  limits  the  ability  of  the  DC  to  properly  “observe”  the  perturbing  effect 
(in  fact,  additional  testing  proves  a  solution  for  the  drag  term  is  not  necessary 
at  this  altitude  for  most  applications;  it  must  be  noted,  however,  that  in  some 
analytic  theories,  the  drag  solve-for  can  be  used  as  a  “free  parameter”  in  the 
DC  to  absorb  errors  stemming  from  truncated  force  models  for  perturbations 
other  than  atmospheric  drag).  The  Cowell^  case  indicates  the  combined  effects 
on  drag  and  SRP  are  less  than  4  meters  for  this  particular  case.  As  mentioned 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 

RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

BL 

24.620212 

5 

185.25 

2.05 

1.03 

.21 

BL  (Npq) 

23.177123 

5 

524.85 

1.92 

0.98 

.20 

CowelGi 

0.6683(E-9) 

5 

0.43849(E-3) 

4.61 

3.11 

.62 

Cowell^ 

0.2995392 

5 

3.7135 

3.74 

2.59 

.52 

CowelG 

7.9674481 

5 

250.17 

3.98 

2.69 

.54 

Cowell;) 

7.8043229 

6 

105.07 

4.62 

3.32 

.55 

SGP4 

23.568135 

3 

182.27 

2.28 

0.95 

.32 

SGP4  (drag) 

19.894640 

7 

358.55 

2.17 

1.41 

.20 

SGP 

556.90745 

3 

2994.5 

1.85 

1.02 

.34 

SGP  (drag) 

773.92103 

3 

6828.3 

2.02 

1.10 

.37 

DSST  (opt) 

19.431369 

4 

31.769 

2.23 

1.41 

.35 

Table  10:  TOPEX  Orbit  115  Minutes  Fit  and  115  Minutes  Predict 


previously,  the  high  altitude  circular  case  represents  an  orbital  regime  in  which 
the  best  performance  from  the  analytic  theories  can  be  expected;  however,  the 
errors  from  the  analytic  theories  are  still  6  times  as  large  as  from  optimized 
DSST  (with  optimized  DSST  running  at  speeds  comparable  to  SGP4).  Fur¬ 
thermore,  the  use  of  drag  solve-for  terms  in  the  analytic  theories  did  not  help  to 
absorb  errors  from  truncations  in  perturbation  models  other  than  atmospheric 
drag. 

MOLNIYA  ORBIT: 

Configuration  Notes:  In  this  case,  CowelGi  has  force  models  matching  those 
in  the  truth  (21x21  GEM10B  geopotential,  Jacchia- Roberts  (Schatten)  drag, 
lunar  solar  point  mass  effects,  SRP,  solve-for  drag  and  SRP  parameters  in  DC). 
CowelG  is  identical  to  CowelGi,  with  the  exception  of  a  reduced  geopotential 
(i.e. ,  4x4).  Cowell;)  differs  from  CowelG  only  in  that  the  exact  drag  and  SRP 
terms  from  the  truth  ephemeris  have  been  hardwired  into  the  run  (they  are  not 
being  solved  for  by  the  DC).  The  annotation  of  Npq  or  “drag”  after  a  theory 
indicates  DC  solve-for  parameters  used  in  the  run. 

In  the  Molniya  case,  geopotential  terms  beyond  the  4x4  configuration  of 
Cowell;)  (up  to  the  21x21  limit  in  the  truth  ephemeris)  introduce  almost  65  m 
worth  of  “pure”  error  over  the  twelve  hour  predict.  The  majority  of  this  error 
can  be  attributed  to  neglected  resonance  terms  (the  Molniya  orbit  is  resonant 
at  the  even  orders;  the  Cowell;)  DC  captured  only  the  2nd  and  4th  resonant 
orders).  If  these  results  are  compared  to  those  of  CowelC,  it  is  clear  this  case 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 

RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

BL 

BL  (Npq) 
CowelGi 

N.A. 

N.A. 

0.4356(E-8) 

6 

0.40299(E-4) 

16.05 

11.94 

1.99 

Cowell 

0.119446(E-2) 

28 

577.20 

43.15 

38.81 

1.39 

Cowell;) 

0.111818(E-2) 

29 

64.329 

39.45 

39.34 

1.36 

SGP4 

1906.6052 

15 

6988.8 

11.19 

9.80 

0.65 

SGP4  (drag) 
SGP 

div.  DC 

1942.0303 

14 

6392.6 

9.18 

8.21 

.59 

SGP  (drag) 

1721.1822 

18 

105170. 

11.74 

10.35 

.58 

DSST  (opt) 

34.803726 

8 

174.37 

11.87 

10.87 

1.36 

Table  11:  Molniya  Orbit  12  Hour  Fit  and  12  Hour  Predict 


also  exhibits  observability  problems.  Specifically,  a  poor  solution  for  the  drag 
parameter  was  obtained  by  the  Cowell  DC.  The  force  models  for  BL  proved 
inadequate  for  this  orbit  in  that  the  satellite  impacted  the  earth  during  the  DC 
processing.  Therefore,  neither  accuracy  nor  timing  metrics  could  be  derived  (in 
fairness,  the  BL  theory  was  neither  developed  nor  intended  for  use  with  this 
orbit  type).  Even  though  SGP  may  not  be  practical  for  this  orbit  type,  its 
performance  was  still  analyzed  (very  poor  performance);  however,  since  SGP4’s 
geopotential  model  had  been  modified  to  include  some  12  hour  resonance  terms 
(a  subset  of  terms  from  the  second  and  fourth  orders),  it  was  tested  (see  Hujsak 
and  Hoots  12).  If  an  attempt  was  made  to  solve  for  the  drag  term,  the  DC 
diverged.  SGP4  results  for  which  the  drag  term  is  not  solved-for  are  extremely 
poor.  On  the  contrary,  the  optimized  DSST  performed  extremely  well  in  terms 
of  accuracy.  In  order  to  provide  a  timing  assessment  that  was  more  realistic, 
further  testing  was  undertaken  to  again  fit  a  Harris-Priester  drag  model  to  the 
Jacchia- Roberts  truth: 

These  results  provide  a  more  realistic  assessment  of  the  blend  of  accuracy 
and  computational  efficiency  which  can  be  obtained  from  optimized  DSST  for 
these  orbits. 

It  is  worth  stating  the  Molniya  represents  the  most  challenging  orbit  analyzed 
in  this  study.  This  is  due  to  the  wide  range  of  perturbing  effects  experienced 
by  this  orbit.  Molniya  orbits  are  characterized  by  low  perigee  heights,  which 
introduce  substantial  atmospheric  and  geopotential  effects  (not  to  mention  that 
most  Molniya  orbits  maintain  repeat  groundtrack  constructs,  which  significantly 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 
RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

CowelGi 

2.78328397(E-4) 

6 

0.12155 

14.90 

12.37 

2.06 

CowelG 

1.19451595(E-3) 

28 

583.21 

41.20 

38.53 

1.38 

DSST  (opt) 

34.804297 

10 

174.64 

13.72 

12.36 

1.24 

Table  12:  Molniya  Orbit  (Real  World)  12  Hour  Fit  and  12  Hour  Predict 


contribute  to  resonance);  in  addition,  the  apogee  heights  of  these  orbits  lead  to 
significant  third  body  and  solar  radiation  pressure  effects.  This  wide  range 
of  perturbative  effects  must  be  considered  in  relation  to  the  high  quality  of 
the  observational  data  available  for  these  orbits.  For  these  reasons,  only  a 
perturbation  theory  with  high  quality  force  models  should  be  considered  to 
determine  the  orbits  of  these  objects.  Further  results  for  this  orbit  class  can  be 
found  in  Fonte21. 

GEOSYNCHRONOUS  ORBIT: 

Configuration  Notes:  In  this  case,  CowelGi  has  force  models  matching  those 
in  the  truth  (21x21  GEM10B  geopotential,  lunar  solar  point  mass  effects,  SRP, 
solve- for  SRP  parameter  in  DC).  CowelG  is  identical  to  CowelGi,  with  the 
exception  of  a  reduced  geopotential  (i.e. ,  4x4).  The  annotation  of  Npq  or  “drag” 
after  a  theory  indicates  DC  solve-for  parameters  used  in  the  run.  As  with  the 
other  cases  (except  decay),  DSST  uses  a  two-card  element  set  guess;  DSST* 
uses  the  osculating  elements  from  which  the  truth  trajectory  was  generated 
(resulting  from  a  Cowell  fit  to  an  SGP4  ephemeris  derived  from  the  two  card 
element  set). 

For  the  geosynchronous  orbit,  sub-meter  differences  arise  between  the  Cow¬ 
ell  21x21  and  4x4  geopotential  configurations.  The  impact  of  the  tesseral  terms 
(mainly  resonance)  can  be  clearly  seen  by  comparing  the  CowelG  results  to 
those  of  the  analytic  theories.  Mathematically,  the  4x4  Cowell  configuration 
completely  captures  resonant  terms  through  the  fourth  order;  the  analytic  the¬ 
ories  (with  the  exception  of  SGP4)  do  not  model  tesseral  terms  (SGP4  models 
a  subset  of  the  tesseral  terms  at  the  second  and  fourth  order  (see  Introduction). 
This  lack  of  tesseral  resonance  modelling,  as  well  as  crude  third  body  models, 
limit  the  analytic  theories  to  kilometer  level  accuracy.  The  optimized  DSST 
DC,  which  contains  a  4x4  geopotential  configuration  like  CowelC,  provides  an 
approximate  15  m  predict  error.  This  small  predict  error  can  be  attributed  to 
truncations  in  the  third  body  short  periodic  model  and  the  neglect  of  solar  ra¬ 
diation  pressure  short  periodics  (truncations  made  to  enhance  DSST’s  speed). 
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Theory 

DC 

RMS  (m) 

ITER 

Predict 

RMS  (m) 

CPU 

(sec) 

DC 

CPU 

CPU 

per 

ITER 

BL 

680.91647 

4 

2152.8 

6.31 

4.96 

1.24 

BL  (Npq) 

682.19405 

4 

19531.0 

7.07 

4.73 

1.18 

CowelGi 

0.2114(E-8) 

5 

0.63256(E-3) 

21.18 

16.20 

3.24 

CowelG 

0.231096(E-2) 

5 

0.2085 

16.18 

12.92 

2.58 

SGP4 

909.57138 

3 

1245.9 

8.24 

7.07 

2.36 

SGP4  (drag) 

673.00834 

3 

9427.0 

8.63 

6.77 

2.26 

SGP 

701.47090 

3 

2020.7 

8.40 

6.01 

2.00 

SGP  (drag) 

701.63614 

3 

2188.7 

7.76 

5.86 

1.95 

DSST*  (opt) 

3.733008 

7 

11.219 

11.45 

9.36 

1.34 

DSST  (opt) 

4.6368485 

9 

10.778 

13.09 

11.15 

1.24 

Table  13:  Geosynchronous  Orbit  1  Day  Fit  and  1  Day  Predict 


The  optimized  DSST  configuration  also  contains  a  weak  time  dependent  for¬ 
mulation,  which  accounts  for  the  movement  of  third  bodies  over  the  averaging 
interval.  A  time  independent  formulation  assumes  the  third  bodies  do  not  move 
over  the  course  of  an  averaging  interval  (typically  one  orbital  period,  which 
for  the  geosynchronous  case  equals  one  day)  (see  Green20).  Obviously,  this  as¬ 
sumption  is  reasonable  for  low  altitude  objects;  however,  this  assumption  breaks 
down  for  objects  at  much  higher  altitudes.  Previous  studies  indicate  not  mod¬ 
elling  weak  time  dependence  adds  about  a  200  m  predict  error  over  the  course 
of  three  days  and  increased  numbers  of  iterations  for  DC  convergence  (see  Fonte 
and  Sabol16). 

In  Table  13,  results  for  two  separate  DSST  DCs  are  listed.  An  inspection  of 
the  number  of  iterations  for  the  DSST  case  emphasizes  the  two-card  element  set 
used  by  the  Space  Command  theories  is  based  upon  a  different  set  of  “mean” 
elements  than  those  used  by  the  Draper  Semianalytic  Satellite  Theory.  DSST* 
uses  initial  conditions  which  are  slightly  more  appropriate  than  the  two  card 
element  set;  however,  Draper  “mean”  elements  still  would  provide  the  best 
initial  conditions  for  the  DC. 
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In  general,  the  orbit  propagators  can  be  ranked  in  terms  of  accuracy: 


Theory 

Decay 

HST 

TOPEX 

Molniya 

Geosynchronous 

Average 

CowelGi 

1 

1 

1 

1 

1 

1 

DSST  (opt) 

2 

2 

2 

2 

3 

2.2 

Cowell 

3 

3 

3 

3 

2 

2.8 

SGP4 

4 

4 

4 

5 

4 

4.2 

BL 

5 

5 

5 

6 

6 

5.4 

SGP 

6 

6 

6 

4 

5 

5.4 

Table  14:  Accuracy  Rankings 

In  consideration  of  these  accuracy  rankings,  it  should  be  noted  the  ana¬ 
lytic  theories  performed  very  poorly  for  the  Molniya  and  geosynchronous  cases; 
predict  accuracies  were  on  the  order  of  kilometers. 


Theory 

Decay 

HST 

TOPEX 

Molniya 

Geosynchronous 

CowelGi  vs  DSST  (opt) 

32.67% 

51.61% 

51.63% 

7.92% 

45.94% 

DSST  (opt)  was 

faster 

faster 

faster 

faster 

faster 

SGP4  vs  DSST  (opt) 

98.05% 

26.05% 

2.19% 

22.61% 

38.96% 

DSST  (opt)  was  ^ 

slower 

slower 

faster 

slower 

slower 

Table  15:  Total  CPU  Timing  Comparisons 

In  consideration  of  the  timing  metrics,  remember  the  optimized  DSST  results 
are  for  particular  configurations  of  DSST  developed  by  Fonte  and  Sabol16.  Their 
rule  of  thumb  in  the  development  of  these  configurations  was  to  keep  DSST  5- 
10  times  more  accurate  than  SGP4.  If  speed  is  truly  an  issue,  these  DSST 
configurations  can  be  further  tailored  to  decrease  DSST’s  run  time;  the  user 
is  free  to  configure  the  theory  at  run  time  to  meet  personal  speed  vs  accuracy 
requirements. 

As  a  final  note,  it  should  be  mentioned  the  timing  statistics  aren’t  exactly 
fair.  As  was  discussed  previously,  a  two-card  element  set  guess  was  used  as  initial 
conditions  for  the  various  differential  correction  runs.  SGP4  and  SGP  have  been 
tailored  to  use  the  “mean”  elements  provided  on  the  two-card  element  set;  DSST 
and  Cowell  DCs,  which  can  use  the  two  card  element  sets  as  initial  conditions, 
would  run  more  efficiently  with  their  own  tailored  initial  conditions  (i.e. ,  Draper 
“mean”  elements  for  DSST  and  osculating  elements  for  Cowell).  In  addition, 

^These  time  metrics  are  based  on  the  short  fit  span  chosen.  Other  studies16  indicate  that 
opt  DSST  performs  more  favorably  with  longer  fit  spans. 
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limited  mismatching  of  force  models  was  done  in  this  analysis  for  Cowell  1  (with 
the  exception  of  using  the  Harris- Priester  model  for  the  drag  perturbed  cases);  in 
essence,  the  Cowell2i  DC  were  simply  attempting  to  recover  appropriate  initial 
conditions.  Operational  scenarios  would  be  less  optimistic. 

Conclusion 


This  paper  evaluated  the  performance  of  various  orbit  propagation  theo¬ 
ries  for  artificial  earth  satellites  in  different  orbital  regimes.  Specifically,  R&D 
GTDS’s  Cowell  (numerical),  DSST  (semianalytical),  SGP,  SGP4,  and  BL  (ana¬ 
lytic)  orbit  propagators  were  compared  for  decaying  circular  (~  200  km  perigee 
height),  low  altitude  circular  (590  km  perigee  height),  high  altitude  circular 
(1340  km  perigee  height),  Molniya,  and  geosynchronous  orbits.  (R&D  GTDS 
was  chosen  because  it  has  all  these  theories  available  in  it.)  Computation  time 
and  RMS  errors  were  used  as  comparison  metrics  on  a  SUN  Sparc  10. 

It  should  be  noted  this  study  was  theoretical  in  nature  (i.e. ,  attempting 
to  understand  the  limitation  of  various  orbit  propagation  theories  with  “GPS- 
like” ,  equally  spaced  observational  data).  Without  question,  special  perturba¬ 
tion  techniques  with  rigorous  force  models  proved  the  most  accurate;  however, 
this  accuracy  comes  at  the  cost  of  computational  efficiency  (an  issue  whose  im¬ 
portance  decreases  as  computing  horsepower  increases  and  becomes  commonly 
available).  As  expected,  the  analytic  theories  (general  perturbation  techniques) 
performed  poorly  in  terms  of  accuracy  (due  to  the  severely  truncated  force  mod¬ 
els),  but  were  efficient.  DSST,  which  represents  a  hybrid  of  special  and  general 
perturbation  techniques,  provided  accuracies  approaching  those  of  special  per¬ 
turbation  techniques  at  speeds  comparable  to  the  analytic  theories. 

Clearly,  improving  the  orbit  determination  process  also  includes  equal  con¬ 
sideration  of  the  observational  data  (to  include  station  coordinates,  biases,  stan¬ 
dard  deviations,  solar  and  geomagnetic  activity,  refraction,  coordinate  frames, 
and  timing  issues).  In  addition,  propagating  an  element  set  in  a  manner  con¬ 
sistent  with  which  it  was  determined  is  also  of  paramount  importance.  With 
the  widespread  technological  advances  of  computing  platforms,  orbit  propaga¬ 
tion  theories,  and  timing  mechanisms,  as  well  as  strong  observational  data  (laser, 
GPS,  radar,  transponder,  etc),  the  time  is  ripe  to  re-analyze  current  operational 
orbit  determination  practices.  Specifically,  an  analysis  should  be  undertaken  to 
trade  the  long  term  cost  implication  of  upgrading  current  orbit  determination 
techniques  to  modern  hardware,  software,  and  astrodynamics  capabilities.  The 
authors  recommend  that  this  study  be  performed  by  a  wide  range  of  organiza¬ 
tions  to  tap  a  broad  source  of  knowledge  and  expertise. 
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